#include "StdAfx.h"
#include "CrossbarNetwork.h"
#include "NetworkBase.h"
#include "StreamIO.h"
#include "PNGData.h"

namespace SpikingNeuronLib
{
	// Public Members -------------------------------------------------------------------------------------------

	CrossbarNetwork::CrossbarNetwork()
	{
		// defaults:
		//int Ne = 800;			// excitatory neurons
		//int Ni = 200;			// inhibitory neurons
		//int M = 100;			// the number of synapses per neuron
		//int D = 20;			// maximal axonal conduction delay

		unsigned int seed = (unsigned int) DateTime::Now.ToFileTimeUtc();
		this->Init(800, 200, 20, seed);
	}

	CrossbarNetwork::CrossbarNetwork(int numExcitatoryNeurons, int numInhibitoryNeurons, int maxDelay)
	{
		unsigned int seed = (unsigned int) DateTime::Now.ToFileTimeUtc();
		this->Init(numExcitatoryNeurons, numInhibitoryNeurons, maxDelay, seed);
	}

	CrossbarNetwork::CrossbarNetwork(int numExcitatoryNeurons, int numInhibitoryNeurons, int maxDelay, unsigned int seed)
	{
		this->Init(numExcitatoryNeurons, numInhibitoryNeurons, maxDelay, seed);
	}

	CrossbarNetwork::~CrossbarNetwork()
	{
	}

	CrossbarNetwork^ CrossbarNetwork::Clone()
	{
		using namespace System::IO;
		using namespace System::Text;

		// save the network state of the current network to a buffer
		MemoryStream^ buffer = gcnew MemoryStream();
		StreamWriter^ sw = gcnew StreamWriter(buffer, Encoding::UTF8, 4096, true);

		try
		{
			this->SaveNetworkState(sw);
			sw->Close();
		}
		finally
		{
			delete sw;
		}

		// load the network state into a new network
		CrossbarNetwork^ network = gcnew CrossbarNetwork(_numExcitatoryNeurons, _numInhibitoryNeurons, _maxDelay);
        buffer->Seek(0L, SeekOrigin::Begin);
        StreamReader^ sr = gcnew StreamReader(buffer, Encoding::UTF8);

		try
		{
			network->LoadNetworkState(sr);
			sr->Close();
		}
		finally
		{
			delete sr;
		}

		return network;
	}

	void CrossbarNetwork::CopyTo(CrossbarNetwork^ network)
	{
		using namespace System::IO;
		using namespace System::Text;

		// save the network state of the current network to a buffer
		MemoryStream^ buffer = gcnew MemoryStream();
		StreamWriter^ sw = gcnew StreamWriter(buffer, Encoding::UTF8, 4096, true);

		try
		{
			this->SaveNetworkState(sw);
			sw->Close();
		}
		finally
		{
			delete sw;
		}

		// load the network state into the other network
        buffer->Seek(0L, SeekOrigin::Begin);
        StreamReader^ sr = gcnew StreamReader(buffer, Encoding::UTF8);

		try
		{
			network->LoadNetworkState(sr);
			sr->Close();
		}
		finally
		{
			delete sr;
		}
	}

	void CrossbarNetwork::SaveFiringsData(String^ filePath, bool includeOverflow)
	{
		using namespace System::IO;

		StreamWriter^ sw = File::CreateText(filePath);

		try
		{
			for (int i = 0; i < _N_firings; i++)
			{
				// include all lines if includeOverflow is set
				// otherwise include the header [-20 0] and all lines with positive firing time
				if (includeOverflow || _firings[i].Time == -20 || _firings[i].Time >= 0)
				{
					sw->WriteLine("{0,3} {1,3}", _firings[i].Time, _firings[i].NeuronIndex);
				}
			}

			sw->Close();
		}
		finally
		{
			delete sw;
		}
	}

	void CrossbarNetwork::LoadFiringsData(String^ filePath)
	{
		using namespace System::IO;

		StreamReader^ sr = File::OpenText(filePath);

		try
		{
			int index = 0;
			while ( sr->Peek() >= 0 )
			{
				for (int j = 0; j < 2; j++)
				{
					Nullable<int> item = StreamIO::ScanStream<int>(sr);
					if (item.HasValue)
					{
						if (j == 0)
							_firings[index].Time = item.Value;
						else
							_firings[index].NeuronIndex = item.Value;
					}
					else if (j == 0)
					{
						index--;
						break;	// empty last line
					}
					else
					{
						String^ message = String::Format("Error in line length. Expecting {0} items on this line.", 2);
						throw gcnew ArgumentException(message);
					}
				}
				index++;
				if (index >= _N_firings_max)
				{
					ResizeFiringStack(_N_firings_max + 10000);
				}
			}

			_N_firings = index;	// next available index
		}
		finally
		{
			delete sr;
		}
	}

	void CrossbarNetwork::SaveConnectivity(String^ filePath)
	{
		using namespace System::IO;

		StreamWriter^ sw = File::CreateText(filePath);

		try
		{
			sw->WriteLine(_connections->Length - 1);
			for (int offset = 1; offset < _connections->Length; offset++)
			{
				Connection^ connection = _connections[offset];
				int excitatory = connection->PreNeuron < _numExcitatoryNeurons;
				sw->WriteLine("{0} {1} {2} {3}", connection->PreNeuron, connection->PostNeuron, connection->Delay, excitatory);
			}

			sw->Close();
		}
		finally
		{
			delete sw;
		}
	}

	void CrossbarNetwork::LoadConnectivity(String^ filePath)
	{
		using namespace System::IO;

		StreamReader^ sr = File::OpenText(filePath);

		try
		{
			int numConnections = StreamIO::ScanStreamForType<int>(sr);
			List<Connection^>^ connections = gcnew List<Connection^>(numConnections);

			for (int offset = 0; offset < numConnections; offset++)
			{
				int preNeuron = StreamIO::ScanStreamForType<int>(sr);
				int postNeuron = StreamIO::ScanStreamForType<int>(sr);
				int delay = StreamIO::ScanStreamForType<int>(sr);
				int isExcitatory = StreamIO::ScanStreamForType<int>(sr);

				connections->Add(gcnew Connection(preNeuron, postNeuron, delay, (isExcitatory > 0)? true : false));
			}

			this->PrepareConnections(connections);
		}
		finally
		{
			delete sr;
		}
	}

	void CrossbarNetwork::SaveNetworkState(StreamWriter^ sw)
	{
		sw->WriteLine(String::Format("V{0,2:F1}", NetworkStateFileVersion));
		sw->WriteLine(_connections->Length - 1);
		for (int offset = 1; offset < _connections->Length; offset++)
		{
			// save the connection, including current synaptic weight and derivative
			Connection^ connection = _connections[offset];
			sw->WriteLine("{0} {1} {2} {3,5:F3} {4,6:F5}", connection->PreNeuron, connection->PostNeuron, connection->Delay, connection->Weight, connection->Derivative);
		}

		// find the earliest spike timestamp within MaxSpikeInteractionDistance of the current timestamp
		long long earliestTimestamp = Int64::MaxValue;
		for (int i = 0; i < _totalNeurons; i++)
		{
			if ((_lastSpikeTime[i] < earliestTimestamp) && (_timestamp - _lastSpikeTime[i] < MaxSpikeInteractionDistance))
				earliestTimestamp = _lastSpikeTime[i];
		}

		for (int i = 0; i < _totalNeurons; i++)
		{
			sw->WriteLine("{0,5:F3} {1,5:F3}", _v[i], _u[i]);
			int spiketime = (int)(_lastSpikeTime[i] - earliestTimestamp);	// values earlier then the earliest value will be negative (these are ignored)
			sw->WriteLine("{0,4:F0}", spiketime);
			sw->WriteLine("{0,9:F8}", _neuronDerivative[i]);
			sw->WriteLine("{0,9:F8}", _depressionAmplitudeScalePrevious[i]);
		}

		sw->Write(" {0} ", _N_firings);
		for (int i = 0; i < _N_firings; i++)
		{
			String^ eventDescription = String::Format("{0} {1} {2} ", _firings[i].Time, _firings[i].NeuronIndex, Enum::GetName(EventLabel::typeid, _firings[i].Label));
			sw->Write(eventDescription);
		}
	}

	void CrossbarNetwork::SaveNetworkState(String^ filePath)
	{
		using namespace System::IO;

		StreamWriter^ sw = File::CreateText(filePath);

		try
		{
			SaveNetworkState(sw);
			sw->Close();
		}
		finally
		{
			delete sw;
		}
	}

	void CrossbarNetwork::LoadNetworkState(StreamReader^ sr)
	{
		double version = StreamIO::CheckVersion(sr);
		int numConnections = StreamIO::ScanStreamForType<int>(sr);
		List<Connection^>^ connections = gcnew List<Connection^>(numConnections);

		// MG TODO load (and also save) the hard and soft weight limits (defaulting to 0 for min, and 10 for max soft and max hard)
		for (int offset = 0; offset < numConnections; offset++)
		{
			int preNeuron = StreamIO::ScanStreamForType<int>(sr);
			int postNeuron = StreamIO::ScanStreamForType<int>(sr);
			int delay = StreamIO::ScanStreamForType<int>(sr);
			double weight = StreamIO::ScanStreamForType<double>(sr);
			double derivative = StreamIO::ScanStreamForType<double>(sr);
			if (preNeuron >= _numExcitatoryNeurons)
			{
				// some older network state files have non-zero derivatives
				// for inhibitory connections due to a bug (now fixed)
				derivative = 0.0;
			}

			connections->Add(gcnew Connection(preNeuron, postNeuron, delay, weight, derivative));
		}

		this->PrepareConnections(connections);

		if (version == NetworkStateFileVersion)
		{
			for (int i = 0; i < _totalNeurons; i++)
			{
				_v[i] = StreamIO::ScanStreamForType<double>(sr);
				_u[i] = StreamIO::ScanStreamForType<double>(sr);

				_lastSpikeTime[i] = (long)StreamIO::ScanStreamForType<int>(sr);
				_neuronDerivative[i] = StreamIO::ScanStreamForType<double>(sr);
				_depressionAmplitudeScalePrevious[i] = StreamIO::ScanStreamForType<double>(sr);
			}
		}
		else
		{
			// reading an old network state format
			// Main differences:
			//	1. Uses the old _LTP[] and _LTD[] plasticity values rather than last spike time
			//	2. Uses boolean flag to denote background event rather than an event label
			// Reconstruct last firing time using the file LTD values (lookup offset into exponential curve)
			// Note this requires that the timestamp get set forward so that last firing times are positive
			_timestamp = MaxSpikeInteractionDistance;
			for (int i = 0; i < _totalNeurons; i++)
			{
				_v[i] = StreamIO::ScanStreamForType<double>(sr);
				_u[i] = StreamIO::ScanStreamForType<double>(sr);

				// approximate the last spike time by comparing the current ltd value with the exponential
				double ltd = StreamIO::ScanStreamForType<double>(sr);
				int index = FindIndex(_LTDLookup, ltd);		// index range 0 to MaxSpikeInteractionDistance (or -1 otherwise)
				_lastSpikeTime[i] = (index > 0) ? _timestamp - index : -1;

				// assume no metaplasticity for old network state format
				_neuronDerivative[i] = 0.0;
				_depressionAmplitudeScalePrevious[i] = 0.0;

				for (int j = 0; j < _maxDelay + 1; j++)
				{
					// consume the values
					double ltp = StreamIO::ScanStreamForType<double>(sr);
				}
			}
		}

		_N_firings = StreamIO::ScanStreamForType<int>(sr);
		if (_N_firings >= _N_firings_max)
		{
			ResizeFiringStack(_N_firings + 10000);
		}

		for (int i = 0; i < _N_firings; i++)
		{
			int time = StreamIO::ScanStreamForType<int>(sr);
			int neuronIndex = StreamIO::ScanStreamForType<int>(sr);
			// Versions:
			//	All versions > 0 use the new FiringEvent with EventLabel (previously a bool)
			//	V1.0 uses new plasticity implementation to support metaplasticity
			//	see NetworkStateFileVersion
			EventLabel label;
			if (version > 0)
			{
				label = StreamIO::ScanStreamForType<EventLabel>(sr);
			}
			else
			{
				bool isBackgroundEvent = StreamIO::ScanStreamForType<bool>(sr);
				label = (isBackgroundEvent) ? EventLabel::Background : EventLabel::Unspecified;
			}
			_firings[i] = FiringEvent(time, neuronIndex, label);
		}	
	}

	void CrossbarNetwork::LoadNetworkState(String^ filePath)
	{
		using namespace System::IO;

		StreamReader^ sr = File::OpenText(filePath);

		try
		{
			LoadNetworkState(sr);
		}
		finally
		{
			delete sr;
		}
	}

	array<FiringEvent>^ CrossbarNetwork::GetFiringsData(bool includeOverflow)
	{
		List<FiringEvent>^ events = gcnew List<FiringEvent>(_N_firings);

		// skip over the firing data with times less than 1
		// (unless includeOverflow is set)
		// as this has been copied over from the previous second
		for (int i = 0; i < _N_firings; i++)
		{
			if (includeOverflow || _firings[i].Time >= 0)
			{
				events->Add(_firings[i]);
			}
		}

		// optionally set the flag on each firing event as foreground, background or internal
		// foreground and background events are externally stimulated
		List<FiringEvent>^ result;
		if (_labelEvents)
		{
			List<FiringEvent>^ labeledEvents = gcnew List<FiringEvent>(events->Count);

			// create a dictionary of FiringEvents keyed by millisecond
			Dictionary<int, List<FiringEvent>^>^ firingEventDictionary = gcnew Dictionary<int, List<FiringEvent>^>();
			for each (FiringEvent firingEvent in events)
			{
				List<FiringEvent>^ currentEvents;
				int time = firingEvent.Time;

				if (firingEventDictionary->ContainsKey(time))
				{
					currentEvents = firingEventDictionary[time];
				}
				else
				{
					currentEvents = gcnew List<FiringEvent>();
					firingEventDictionary->Add(time, currentEvents);
				}

				currentEvents->Add(firingEvent);
			}

			// 8 msec windows that label background and foreground firing events
			// i.e. they associate a firing event with background or foreground stimulation
			array<unsigned char>^ backgroundAssociationWindow = gcnew array<unsigned char>(_totalNeurons);
			array<unsigned char>^ foregroundAssociationWindow = gcnew array<unsigned char>(_totalNeurons);

			// run through the last second
			// if an external (foreground or background) event was fired then start a countdown on each fired neuron
			// if a FiringEvent is generated in the following milliseconds for those neurons
			// then they can be labelled as foreground or background
			// the current millisecond is skipped by labeling FIRST and then starting the countdown on external events
			for (int t = 0; t < 1000; t++)		// simulation of 1 sec
			{
				// any events in the current millisecond?
				if (firingEventDictionary->ContainsKey(t))
				{
					List<FiringEvent>^ firingEvents = firingEventDictionary[t];

					for each (FiringEvent firingEvent in firingEvents)
					{
						int neuronIndex = firingEvent.NeuronIndex;

						// select the label: foreground, background or internal
						EventLabel label = EventLabel::Unspecified;
						if (foregroundAssociationWindow[neuronIndex])
						{
							foregroundAssociationWindow[neuronIndex] = 0;
							label = EventLabel::Foreground;
						}
						else if (backgroundAssociationWindow[neuronIndex])
						{
							backgroundAssociationWindow[neuronIndex] = 0;
							label = EventLabel::Background;
						}
						else
						{
							// default: an internal event
							label = EventLabel::Internal;
						}

						FiringEvent labeledEvent = FiringEvent(firingEvent.Time, firingEvent.NeuronIndex, label);
						labeledEvents->Add(labeledEvent);
					}
				}

				// decrement the association windows
				for (int i = 0; i < _totalNeurons; i++)
				{
					if (backgroundAssociationWindow[i])
						backgroundAssociationWindow[i]--;
					if (foregroundAssociationWindow[i])
						foregroundAssociationWindow[i]--;
				}

				if (_stimulusFiringData && _stimulusFiringData[t])
				{
					array<int>^ selectedNeurons = _stimulusFiringData[t];
					for each (int neuron in selectedNeurons)
					{
						// start a countdown for association between this stimulation
						// and any subsequent firing events (default = 8 msec)
						if (!foregroundAssociationWindow[neuron])
							foregroundAssociationWindow[neuron] = DefaultAssociationWindowSize;
					}
				}

				if (_backgroundFiringData && _backgroundFiringData[t])
				{
					array<int>^ selectedNeurons = _backgroundFiringData[t];
					for each (int neuron in selectedNeurons)
					{
						// start a countdown for association between this stimulation
						// and any subsequent firing events (default = 8 msec)
						if (!backgroundAssociationWindow[neuron])
							backgroundAssociationWindow[neuron] = DefaultAssociationWindowSize;
					}
				}
			}

			result = labeledEvents;
		}
		else
		{
			result = events;
		}

		return result->ToArray();
	}

	// For grid networks, get all of the connections from the specified pre-synaptic neuron that have the specified delay
	array<int>^ CrossbarNetwork::GetConnectionsForNeuronAndDelay(int neuronIndex, int delay)
	{
		// delay is stored as zero-based so delays 1-20 msec map to index range 0-19
		if ((_connectionsWithDelay != nullptr) && (neuronIndex >= 0) && (neuronIndex < _connectionsWithDelay->GetLength(0)) && (delay > 0) && (delay <= _connectionsWithDelay->GetLength(1)))
		{
			return _connectionsWithDelay[neuronIndex, delay - 1];
		}
		else
		{
			return gcnew array<int>(0);
		}
	}

	// Set all connection weights to a random value up to the specified maximum weight
	void CrossbarNetwork::RandomizeWeights(double maxWeight)
	{
		unsigned int seed = (unsigned int)DateTime::Now.ToFileTimeUtc();
		Random^ random = gcnew Random(seed);

		for each (Connection^ connection in _connections)
		{
			if (connection != nullptr)	// first connection slot is a dummy slot
			{
				if (connection->PreNeuron < _numExcitatoryNeurons)
				{
					connection->Weight = maxWeight * random->NextDouble();
				}
				else
				{
					connection->Weight = -(maxWeight * random->NextDouble());
				}
			}
		}

		random = nullptr;
	}

	// Perturb weight values up or down in a random direction and by a random amount within the specified range
	void CrossbarNetwork::PerturbWeights(int perturbationRange)
	{
		unsigned int seed = (unsigned int)DateTime::Now.ToFileTimeUtc();
		Random^ random = gcnew Random(seed);

		for each (Connection^ connection in _connections)
		{
			if (connection != nullptr)	// first connection slot is a dummy slot
			{
				double direction = random->NextDouble() > 0.5 ? -1.0 : 1.0;
				double perturbation = perturbationRange * random->NextDouble();
				double newWeightValue = connection->Weight + (direction * perturbation);

				if (connection->PreNeuron < _numExcitatoryNeurons)
				{
					connection->Weight = newWeightValue >= 0.0 ? newWeightValue : 0.0;
				}
				else
				{
					connection->Weight = newWeightValue <= 0.0 ? newWeightValue : 0.0;
				}
			}
		}

		random = nullptr;
	}

	// Set all excitatory and inhibitory connection weights to the specified values
	void CrossbarNetwork::SetUniformConnectionWeights(double excitatoryWeight, double inhibitoryWeight)
	{
		for each (Connection^ connection in _connections)
		{
			if (connection != nullptr)	// first connection slot is a dummy slot
			{
				if (connection->PreNeuron < _numExcitatoryNeurons)
				{
					connection->Weight = excitatoryWeight;
				}
				else
				{
					connection->Weight = inhibitoryWeight;
				}
			}
		}
	}

	void CrossbarNetwork::AddFixedConnectionList(int synapsesPerNeuron)
	{
		this->AddConnections(this->CreateFixedConnectionList(_numExcitatoryNeurons, _numInhibitoryNeurons, synapsesPerNeuron, _maxDelay, 6.0, -5.0));
	}

	void CrossbarNetwork::AddFixedConnectionList(int synapsesPerNeuron, double defaultExcitatoryWeight, double defaultInhibitoryWeight)
	{
		this->AddConnections(this->CreateFixedConnectionList(_numExcitatoryNeurons, _numInhibitoryNeurons, synapsesPerNeuron, _maxDelay, defaultExcitatoryWeight, defaultInhibitoryWeight));
	}

	void CrossbarNetwork::AddConnections(List<Connection^>^ connections)
	{
		this->PrepareConnections(connections);
	}

	void CrossbarNetwork::AddConnection(Connection^ connection)
	{
		if (_singles == nullptr)
		{
			_singles = gcnew List<Connection^>();
		}

		_singles->Add(connection);

		// update connection numbers
		_numConnections = _singles->Count;					// the total number of connections in the network
		_numExcitatoryConnections = 0;
		_numInhibitoryConnections = 0;
		for each (Connection^ connection in _singles)
		{
			if (connection->PreNeuron < _numExcitatoryNeurons)
			{
				_numExcitatoryConnections++;
			}
			else
			{
				_numInhibitoryConnections++;
			}
		}
	}

	// Private Members -------------------------------------------------------------------------------------------

	void CrossbarNetwork::Init(int numExcitatoryNeurons, int numInhibitoryNeurons, int maxDelay, unsigned int seed)
	{
		_initialized = false;

		Random^ random = gcnew Random(seed);

		_numExcitatoryNeurons = numExcitatoryNeurons;							// excitatory neurons
		_numInhibitoryNeurons = numInhibitoryNeurons;							// inhibitory neurons
		_totalNeurons = _numExcitatoryNeurons + _numInhibitoryNeurons;
		_maxDelay = maxDelay;													// maximal axonal conduction delay

		_N_firings_max = 15000;					// Izhikevich uses 100,000
		_minSynapticWeightSoftLimit = 0.0;		// same as the implicit hard limit for min synaptic weight
		_maxSynapticWeightSoftLimit = C_max;	// soft limit == hard limit for max synaptic weight (by default)
		_maxSynapticWeightHardLimit = C_max;
		_potentiationAmplitudeBaseline = PotentiationAmplitudeDefault;	// default value for STDP rule constant for LTP window = 0.1
		_depressionAmplitudeBaseline = DepressionAmplitudeDefault;		// default value for STDP rule constant for LTD window = 0.12
		_foregroundStimulationLevel = ForegroundStimulationLevelDefault;
		_backgroundStimulationLevel = BackgroundStimulationLevelDefault;

		// Initialise metaplasticity variables
		_potentiationAmplitudeScale = gcnew array<double>(_totalNeurons);
		_depressionAmplitudeScale = gcnew array<double>(_totalNeurons);
		_depressionAmplitudeScalePrevious = gcnew array<double>(_totalNeurons);		// scale values for each neuron in the previous second
		_inertia = MetaplasticityDefaultInertia;
		_upperWeightLimitResistance = MetaplasticityDefaultWeightLimitResistance;
		_lowerWeightLimitResistance = MetaplasticityDefaultWeightLimitResistance;
		_upperWeightLimitPrecision = MetaplasticityDefaultWeightLimitPrecision;
		_lowerWeightLimitPrecision = MetaplasticityDefaultWeightLimitPrecision;

		_lastSpikeTime = gcnew array<long long>(_totalNeurons);
		for (int i = 0; i < _totalNeurons; i++)
			_lastSpikeTime[i] = -1;

		_wsd = gcnew array<double>(_totalNeurons);
		_neuronDerivative = gcnew array<double>(_totalNeurons);
		_connectionsPerNeuron = gcnew array<int>(_totalNeurons);

		// create fixed plasticity lookup tables for exponential decay over 1000 msec
		// MG this is an approximation of an exponential from the original code; replace this approximation with exponential
		_LTPLookup = gcnew array<double>(MaxSpikeInteractionDistance);
		_LTDLookup = gcnew array<double>(MaxSpikeInteractionDistance);
		_LTPLookup[0] = _potentiationAmplitudeBaseline;
		_LTDLookup[0] = _depressionAmplitudeBaseline;
		for (int i = 1; i < 1000; i++)
		{
			_LTPLookup[i] = 0.95 * _LTPLookup[i - 1];			// tau = 20 ms
			_LTDLookup[i] = 0.95 * _LTDLookup[i - 1];			// tau = 20 ms
		}

		_a = gcnew array<double>(_totalNeurons);
		for (int i = 0; i < _numExcitatoryNeurons; i++)
			_a[i] = 0.02;

		for (int i = _numExcitatoryNeurons; i < _totalNeurons; i++)
			_a[i] = 0.1;

		_d = gcnew array<double>(_totalNeurons);
		for (int i = 0; i < _numExcitatoryNeurons; i++)
			_d[i] = 8;

		for (int i = _numExcitatoryNeurons; i < _totalNeurons; i++)
			_d[i] = 2;

		// initial values for v
		_v = gcnew array<double>(_totalNeurons);
		for (int i = 0; i < _totalNeurons; i++)
			_v[i] = -65 + 10 * random->NextDouble();		// -65.0 inclusive to -55.0 exclusive

		// initial values for u
		_u = gcnew array<double>(_totalNeurons);
		for (int i = 0; i < _totalNeurons; i++)
			_u[i] = 0.2 * _v[i];

		// spike timings
		_firings = gcnew array<FiringEvent>(_N_firings_max);
		_firings[0] = FiringEvent(-_maxDelay, 0, EventLabel::Unspecified);
		_N_firings = 1;

		random = nullptr;
	}

	void CrossbarNetwork::PrepareConnections(List<Connection^>^ connections)
	{
		// prepare a 2D array of lists
		// dimension 1 = pre-synaptic neuron, dimension 2 = the axonal delay
		// subtract one from delay (1 to _maxDelay inclusive) to get 0-based dimension 2
		array<List<int>^, 2>^ connectionDelayOffsets = gcnew array<List<int>^, 2>(_totalNeurons, _maxDelay);
		for (int i = 0; i < _totalNeurons; i++)
		{
			for (int j = 0; j < _maxDelay; j++)
			{
				// create a list of connection offsets
				connectionDelayOffsets[i, j] = gcnew List<int>(5);
			}
		}

		// assume that neuron indices have been supplied such that
		//		indices start from zero
		//		indices below _numExcitatoryNeurons are excitatory
		//		indices above _numExcitatoryNeurons are inhibitory
		// a check is made below to ensure that connections issuing from each neural type
		// have appropriate synaptic weights i.e. positive weight for excitatory neurons, and negative for inhibitory
		_connections = gcnew array<Connection^>(connections->Count + 1);		// extra element; first offset is a dummy
		_crossbar = gcnew array<int, 2>(_totalNeurons, _totalNeurons);
		_numConnections = 0;					// the total number of connections in the network
		_numExcitatoryConnections = 0;
		_numInhibitoryConnections = 0;

		// offsets must all be non-zero (on the crossbar, zero means no connection)
		// therefore, start at one
		for (int offset = 1; offset < _connections->Length; offset++)
		{
			Connection^ connection = connections[offset - 1];
			_connections[offset] = connection;

			if (connection->PreNeuron < _numExcitatoryNeurons)
			{
				_numExcitatoryConnections++;
			}
			else
			{
				_numInhibitoryConnections++;
			}

			// Note that soft limit == hard limit for max synaptic weight (by default)
			if (connection->PreNeuron < 0 || connection->PreNeuron >= _totalNeurons
				|| connection->PostNeuron < 0 || connection->PostNeuron >= _totalNeurons
				|| connection->Weight < -_maxSynapticWeightSoftLimit || connection->Weight > _maxSynapticWeightSoftLimit
				|| connection->Delay < 1 || connection->Delay > _maxDelay)
			{
				throw gcnew ArgumentException(String::Format("Connection parameters out of range: {0}", connection->ToString()));
			}

			if (connection->PreNeuron < _numExcitatoryNeurons && connection->Weight < 0.0)
			{
				throw gcnew ArgumentException(String::Format("Excitatory connections must have positive weight: {0}", connection->ToString()));
			}

			if (connection->PreNeuron >= _numExcitatoryNeurons && connection->Weight > 0.0)
			{
				throw gcnew ArgumentException(String::Format("Inhibitory connections must have negative weight: {0}", connection->ToString()));
			}

			List<int>^ offsets = connectionDelayOffsets[connection->PreNeuron, connection->Delay - 1];
			offsets->Add(offset);

			_crossbar[connection->PreNeuron, connection->PostNeuron] = offset;
		}

		_numConnections = connections->Count;

		// transform lists to arrays
		_connectionsWithDelay = gcnew array<array<int>^, 2>(_totalNeurons, _maxDelay);
		for (int i = 0; i < _totalNeurons; i++)
		{
			for (int j = 0; j < _maxDelay; j++)
			{
				_connectionsWithDelay[i, j] = connectionDelayOffsets[i, j]->ToArray();
			}
		}
		connectionDelayOffsets = nullptr;

		// compute the number of synaptic inputs for each postsynaptic neuron
		// needed for metaplasticity calculation
		for (int offset = 1; offset < _connections->Length; offset++)
		{
			Connection^ connection = _connections[offset];

			_connectionsPerNeuron[connection->PostNeuron] += 1;
		}

		_initialized = true;
	}

	List<Connection^>^ CrossbarNetwork::CreateFixedConnectionList(int numExcitatory, int numInhibitory, int synapsesPerNeuron, int maxDelay, double defaultExcitatoryWeight, double defaultInhibitoryWeight)
	{
		// Note that if synapsesPerNeuron is less than maxDelay then connection lengths will not be distributed over all possible delays
		// i.e. if synapsesPerNeuron = 10 and maxDelay = 20 then connections will have delays 1 to 10 (longer delays will not occur)
		// N1000 networks are not affected as synapsesPerNeuron = 100 and maxDelay = 20 and therefore all delay bands are filled
		int delaybandWidth = synapsesPerNeuron / maxDelay;		// width of each delay band (a band of offsets, each with the same delay)
		if (delaybandWidth == 0) delaybandWidth = 1;

		int totalNeurons = numExcitatory + numInhibitory;
		List<Connection^>^ connections = gcnew List<Connection^>(totalNeurons * synapsesPerNeuron);

		Random^ random = gcnew Random((int)DateTime::Now.ToFileTimeUtc());

		array<int>^ postNeurons = gcnew array<int>(synapsesPerNeuron);		// the post-synaptic neurons for the current pre-synaptic neuron
		for (int pre = 0; pre < totalNeurons; pre++)
		{
			double weight;
			int delay = 0;

			for (int j = 0; j < synapsesPerNeuron; j++) 
			{
				int post;
				bool exists;
				do
				{
					exists = false;
					if (pre < numExcitatory)
					{
						// excitatory neurons can connect to any other neuron
						post = random->Next(totalNeurons);		// 0 to totalNeurons-1
					}
					else
					{
						// inhibitory neurons can only connect to excitatory neurons
						post = random->Next(numExcitatory);		// 0 to numExcitatory-1
					}

					// prevent self connections
					if (post == pre)
					{
						exists = true;
						continue;
					}

					// prevent duplicate connections
					for (int k = 0; k < j; k++)
					{
						if (postNeurons[k] == post)
						{
							exists = true;
							break;
						}
					}
				}
				while (exists);

				if (pre < numExcitatory)
				{
					weight = defaultExcitatoryWeight;	//6.0;

					// excitatory connections have axonal delay of 1 to maxDelay msec
					if (j % delaybandWidth == 0)
					{
						delay++;						// move to the next delay band
						if (delay % maxDelay == 1)		// allow wrap on overflow
						{
							delay = 1;
						}
					}
				}
				else
				{
					weight = defaultInhibitoryWeight;	//-5.0;

					delay = 1;		// all inhibitory connections have axonal delay of 1 msec
				}

				postNeurons[j] = post;
				connections->Add(gcnew Connection(pre, post, delay, weight));
			}
		}

		return connections;
	}

	int CrossbarNetwork::RunImpl(int totalSeconds, PatternGenerator^ backgroundPattern, IStimulus^ inputPattern)
	{
		return RunNetwork(0, totalSeconds, backgroundPattern, inputPattern);
	}

	int CrossbarNetwork::RunImpl(int startSeconds, int endSeconds, PatternGenerator^ backgroundPattern, IStimulus^ inputPattern)
	{
		return RunNetwork(startSeconds, endSeconds, backgroundPattern, inputPattern);
	}

	int CrossbarNetwork::RunNetwork(int startSeconds, int endSeconds, PatternGenerator^ backgroundPattern, IStimulus^ inputPattern)
	{
		int result = 1;

		// Spike interaction model is nearest neighbour with a maximum separation of 1000 msec
		// Note that the original Izhikevich implementation (using _LTP[i, t + _maxDelay] and _LTD[i])
		// allowed unlimited separation by copying over current plasticity values across each second
		// However this implementation was incompatible with the new mechanism to support metaplasticity

		// initialise per neuron plasticity values
		for (int i = 0; i < _totalNeurons; i++)
		{
			_potentiationAmplitudeScale[i] = 1.0;
			_depressionAmplitudeScale[i] = 1.0;
		}

		// initialise any connections added one-by-one
		if (_singles != nullptr)
		{
			this->PrepareConnections(_singles);
			_singles = nullptr;
		}

		// check that the stimulus is compatible with the network
		// i.e. each neuron that fires in the stimulus must exist in the network
		if (inputPattern)
		{
			HashSet<int>^ firingNeurons = inputPattern->PatternNeurons;
			for each (int neuronIndex in firingNeurons)
			{
				if (neuronIndex >= this->TotalNeurons)
				{
					_initialized = false;
					throw gcnew System::ArgumentException("The stimulus references neuron indices outside of the range for this network", "inputPattern");
				}
			}
		}

		if (_initialized)
		{
			// input vector
			array<double>^ I = gcnew array<double>(_totalNeurons);
			_timestamp = -1;

			for (int sec = startSeconds; sec < endSeconds; sec++)
			{
				// get the external stimulus firing times for the current second (both input and background)
				_stimulusFiringData = nullptr;
				if (inputPattern)
				{
					_stimulusFiringData = inputPattern->GetNeuronFiringData(sec);
				}
				_backgroundFiringData = nullptr;
				if (backgroundPattern)
				{
					_backgroundFiringData = backgroundPattern->GetNeuronFiringData(sec);
				}

				for (int t = 0; t < 1000; t++)		// simulation of 1 sec
				{
					_timestamp++;

					// input vector initialisation
					for (int i = 0; i < _totalNeurons; i++)
						I[i] = 0;

					for (int i = 0; i < _totalNeurons; i++)
					{
						if (_v[i] >= 30)
						{
							// neuron i has fired and needs resetting
							_v[i] = -65;
							_u[i] += _d[i];

							_lastSpikeTime[i] = _timestamp;

							// check if any pre-synaptic neurons have recently fired
							for (int j = 0; j < _totalNeurons; j++)
							{
								// pre-before-post firing: the post-synaptic event has just occurred
								// for each connection between pre-synaptic neuron j and post-synaptic neuron i ...
								// MG if j is inhibitory (>=800) the derivative is still incremented!!
								if (_crossbar[j, i] > 0)	// connections on row i
								{
									// adjust the derivative of pre-synaptic neurons (excitatory only)
									if ((j < _numExcitatoryNeurons) && (_lastSpikeTime[j] >= 0))
									{
										int offset = _crossbar[j, i];

										long long spikeArrivalTime = _lastSpikeTime[j] + _connections[offset]->Delay;
										long long index = _timestamp - spikeArrivalTime;
										if (index >= 0 && index < MaxSpikeInteractionDistance)
										{
											// the derivative is about to change
											UpdateSynapticIntegral(offset, t);

											double delta = _LTPLookup[(int)index] * _potentiationAmplitudeScale[i];
											_connections[offset]->Derivative += delta;
										}
									}
								}
							}

							// add new firing event
							if (_N_firings >= _N_firings_max)
							{
								ResizeFiringStack(_N_firings_max + 10000);
							}

							_firings[_N_firings++] = FiringEvent(t, i, EventLabel::Unspecified);
						}
					}

					int k = _N_firings - 1;

					// check for arriving spikes
					// check spike arrival for each event in the last frame
					// where a frame is those events on the stack that fired in the last _maxDelay msec
					while (t - _firings[k].Time < _maxDelay)
					{
						int firedTime = _firings[k].Time;
						int firedNeuron = _firings[k].NeuronIndex;		// a pre-synaptic neuron has recently fired

						// post-synaptic spike arrival occurs on these outgoing connections
						// _connectionsWithDelay is 0-based so: 
						//		[0,0] is firedNeuron 0 with delay 1
						//		[5,19] is firedNeuron 5 with delay 20
						array<int>^ connectionOffsets = _connectionsWithDelay[firedNeuron, t - firedTime];
						for each (int offset in connectionOffsets)
						{
							int postNeuron = _connections[offset]->PostNeuron;
							// internal input that contributes to the possible firing of postNeuron in the next msec
							I[postNeuron] += _connections[offset]->Weight;

							if (firedNeuron < _numExcitatoryNeurons)
							{
								// post-before-pre firing
								// a spike from an excitatory presynaptic neuron has arrived
								// assuming that the postsynaptic neuron has recently fired
								// and therefore that the spike did not contribute to postsynaptic firing
								// depress the synapse according to the plasticity value at the post-synaptic neuron

								if (_lastSpikeTime[postNeuron] >= 0)
								{
									long long index = _timestamp - _lastSpikeTime[postNeuron];
									if (index >= 0 && index < MaxSpikeInteractionDistance)
									{
										// check if the postsynaptic firing was in the current second
										// note possible loss of data in truncating lastSpikeTime / 1000 to a double
										// however 'sec' would overflow before lastSpikeTime got that large
										// (max integer for a 64 bit double is approx. 2^52 or 10^15)
										int postFiringSecond = (int) Math::Floor((double)(_lastSpikeTime[postNeuron] / 1000));

										// select the scale
										double depressionAmplitudeScale;
										if (postFiringSecond < sec)
											depressionAmplitudeScale = _depressionAmplitudeScalePrevious[postNeuron];
										else
											depressionAmplitudeScale = _depressionAmplitudeScale[postNeuron];

										// the derivative is about to change
										UpdateSynapticIntegral(offset, t);

										double delta = _LTDLookup[(int)index] * depressionAmplitudeScale;
										_connections[offset]->Derivative -= delta;
									}
								}
							}
						}

						k--;
					}

					// set external input (internal input has already been set)
					if (_stimulusFiringData)
					{
						// set input vector values (corresponding to specific neurons) for the current millisecond
						if (_stimulusFiringData[t])
						{
							array<int>^ selectedNeurons = _stimulusFiringData[t];
							for each (int neuron in selectedNeurons)
							{
								// reset _u[neuron] to default so that the neuron will fire even if it has recently fired
								_u[neuron] = -13.0;

								// add whatever value is required to get _v[neuron] to the _foregroundStimulationLevel
								double current = (0.04 * _v[neuron] + 5) * _v[neuron] + 140 - _u[neuron] + I[neuron];
								if (current < _foregroundStimulationLevel)
									I[neuron] += (_foregroundStimulationLevel - current);
							}
						}
					}

					// optionally include random background firing
					if (_backgroundFiringData)
					{
						if (_backgroundFiringData[t])
						{
							array<int>^ selectedNeurons = _backgroundFiringData[t];
							for each (int neuron in selectedNeurons)
							{
								I[neuron] += _backgroundStimulationLevel;
							}
						}
					}

					for (int i = 0; i < _totalNeurons; i++)
					{
						_v[i] += 0.5 * ((0.04 * _v[i] + 5) * _v[i] + 140 - _u[i] + I[i]);	// for numerical stability
						_v[i] += 0.5 * ((0.04 * _v[i] + 5) * _v[i] + 140 - _u[i] + I[i]);	// time step is 0.5 ms

						_u[i] += _a[i] * (0.2 * _v[i] - _u[i]);
					}

					OneMillisecondTick(sec, t);

				}	// end of loop over 1000 msecs

				// average firing rates for excitatory neurons and inhibitory neurons
				int excitatoryFiringCount = 0;
				int inhibitoryFiringCount = 0;
				for (int i = 1; i < _N_firings; i++)
				{
					if (_firings[i].Time >= 0)
					{
						if (_firings[i].NeuronIndex < _numExcitatoryNeurons)
							excitatoryFiringCount++;
						else
							inhibitoryFiringCount++;
					}
				}
				double frate = (double)excitatoryFiringCount / _numExcitatoryNeurons;
				double ifrate = (double)inhibitoryFiringCount / _numInhibitoryNeurons;

				// proportion of excitatory-excitatory connections that are near maximal synaptic weight (the soft limit if metaplasticity is enabled)
				int nearMaximal = 0;
				for (int i = 1; i <= _numConnections; i++)
				{
					if ((_connections[i]->PreNeuron < _numExcitatoryNeurons) && (_connections[i]->PostNeuron < _numExcitatoryNeurons) && (_connections[i]->Weight > 0.9 * _maxSynapticWeightSoftLimit))
					{
						nearMaximal++;
					}
				}
				double strength = 100 * (double)nearMaximal / _numExcitatoryConnections;

				// update plasticity variables for each neuron
				if (_metaplasticityEnabled)
				{
					// clear Weighted Synaptic Derivative (WSD) accumulator
					for (int i = 0; i < _totalNeurons; i++)
						_wsd[i] = 0.0;

					for (int offset = 1; offset < _connections->Length; offset++)
					{
						// final update of the synaptic integral for the current second
						UpdateSynapticIntegral(offset, 999);

						Connection^ connection = _connections[offset];
						if (connection->PreNeuron < _numExcitatoryNeurons)
						{
							// resist changes in synaptic weight near the weight limits (minimum and maximum synaptic weight)
							// weight is scaled by the closeness of the weight to the weight limits
							// scaling is computed from a 2D surface function that is parameterized by weight and synaptic derivative

							// compress the derivative into the range 0-10
							double compressedDerivative = 0.5 * (connection->Derivative + 10);
							if (compressedDerivative < 0.0)
								compressedDerivative = 0.0;
							else if (compressedDerivative > 10.0)
								compressedDerivative = 10.0;

							// weighting function 2
							double weightedDerivative =
								(_upperWeightLimitResistance * Math::Exp(_upperWeightLimitPrecision * compressedDerivative * (connection->Weight - _minSynapticWeightSoftLimit))) -
								(_lowerWeightLimitResistance * Math::Exp(_lowerWeightLimitPrecision * (10 - compressedDerivative) * (_maxSynapticWeightSoftLimit - connection->Weight)));

							_wsd[connection->PostNeuron] += weightedDerivative;
						}

						connection->Integral = 0.0;
					}

					// compute the neuron derivative
					// a range-limited average of the synaptic inputs for each neuron (averaged over the synaptic inputs for each neuron)
					for (int i = 0; i < _numExcitatoryNeurons; i++)
					{
						// for neurons with synaptic inputs...
						if (_connectionsPerNeuron[i])
						{
							_depressionAmplitudeScalePrevious[i] = _depressionAmplitudeScale[i];

							// update neuron derivative
							_neuronDerivative[i] = Math::Tanh((_inertia * _wsd[i]) / _connectionsPerNeuron[i]);

							// update plasticity amplitude scale (Neuron Derivative has range -1 to 1)
							_potentiationAmplitudeScale[i] = 1.0 - _neuronDerivative[i];		// range: 0 - 2
							_depressionAmplitudeScale[i] = 1.0 + _neuronDerivative[i];			// range: 0 - 2
						}
					}
				}

				OneSecondTickBeforeCleanup(sec, frate, strength, ifrate);

				if (backgroundPattern)
					backgroundPattern->Reset();		// prepare background firing patterns for the next second

				// find the start of the last frame
				int k = _N_firings - 1;
				while (1000 - _firings[k].Time < _maxDelay)
				{
					k--;
				}

				// copy the last frame to the start of the firing stack
				// adjust firing time
				for (int i = 1; i < _N_firings - k; i++)
				{
					_firings[i].Time = _firings[k + i].Time - 1000;
					_firings[i].NeuronIndex = _firings[k + i].NeuronIndex;
				}

				// truncate the firing stack
				_N_firings = _N_firings - k;

				// clip the weights for each excitatory connection
				// note that inhibitory weights are not plastic
				for (int i = 1; i <= _numConnections; i++)	// first element is a dummy
				{
					if (_connections[i]->PreNeuron < _numExcitatoryNeurons)
					{
						_connections[i]->Derivative *= 0.9;				// tau = 250 ms

						if (!_plasticityDisabled)
						{
							_connections[i]->Weight += 0.01 + _connections[i]->Derivative;
							if (_connections[i]->Weight > _maxSynapticWeightHardLimit)
								_connections[i]->Weight = _maxSynapticWeightHardLimit;
							if (_connections[i]->Weight < 0)
								_connections[i]->Weight = 0;
						}
					}
				}

				// SaveNetworkState() can be called from within this callback, or at any other time
				// (except from within the beforeCleanupStatusCallback as the above cleanup code must be called first)
				OneSecondTickAfterCleanup(sec);

			}	// end loop over many seconds

			delete I;
			I = nullptr;

			result = 0;
		}

		return result;
	}

	void CrossbarNetwork::ResizeFiringStack(int newSize)
	{
		// copy the firing stack to a new larger array
		array<FiringEvent>^ temp = gcnew array<FiringEvent>(newSize);
		for (int evt = 0; evt < _N_firings; evt++)
		{
			temp[evt] = _firings[evt];
		}

		_firings = temp;
		_N_firings_max = _firings->Length;
	}

	// Search through an exponential lookup table to find the index
	// of the value closest to the specified value
	// Used to estimate last spike time from the current LTD value
	// Returns -1 if the value is more than the maximimum interaction distance along the curve
	int CrossbarNetwork::FindIndex(array<double>^ data, double value)
	{
		int index = -1;

		for (int i = 0; i < MaxSpikeInteractionDistance && i < data->Length; i++)
		{
			if (value >= data[i])
			{
				index = i;
				break;
			}
		}

		return index;
	}

	// update the synaptic integral for the specified Connection
	// call this BEFORE the derivative is updated
	inline void CrossbarNetwork::UpdateSynapticIntegral(int offset, int t)
	{
		Connection^ connection = _connections[offset];

		// number of timesteps since last delta in the current second
		int timesteps = (t > connection->LastDelta) ?
			// last delta in the current second
			t - connection->LastDelta :
			// last delta occurred in the previous second or earlier
			// just count steps in the current second
			timesteps = t;

		connection->Integral += (connection->Derivative * timesteps);
		connection->LastDelta = t;
	}
}
